use dataset_with_independents, clear

* Collapse

collapse ft_dist polar_dif, by(year)
sum year
gen time = (year-r(min))/(r(max)-r(min))
orthpoly time, generate(time1 time2) degree(2)

* Bimodality

twoway ///
(scatter polar_dif year) ///
(fpfit polar_dif year) ///
, ///
title("Bimodality") ///
xtitle("Year") ///
ytitle("Esteban--Ray polarization index") ///
ylab(, angle(horiz) nogrid) ///
legend(rows(1) order(2 1) label(2 "Estimates") label(1 "Polynomial trend")) ///
scheme(s2mono) ///
graphregion(fcolor(white) lcolor(white)) 

gr export "figureA4.eps", replace
